RTS Model V2 Performance Analysis

TODO:

  • Figure out why there appear to be non-zero IoU scores for predictions that didn’t find any RTS.
  • Possible to remove extra influence of multiple RTS within one tile?
  • Figure out why my IoU scores don’t quite match Yili’s - Check removing non-existant predictions.
  • Add legend for the validation and prediction polygons
  • Shapley analysis (figure out how to import Python models from Google Drive)

Set-Up

Load Libraries

## /home/hrodenhizer/miniconda3/envs/planet_processing_test/lib/python3.8/site-packages/shap/utils/_clustering.py:35: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
##   def _pt_shuffle_rec(i, indexes, index_mask, partition_tree, M, pos):
## /home/hrodenhizer/miniconda3/envs/planet_processing_test/lib/python3.8/site-packages/shap/utils/_clustering.py:54: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
##   def delta_minimization_order(all_masks, max_swap_size=100, num_passes=2):
## /home/hrodenhizer/miniconda3/envs/planet_processing_test/lib/python3.8/site-packages/shap/utils/_clustering.py:63: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
##   def _reverse_window(order, start, length):
## /home/hrodenhizer/miniconda3/envs/planet_processing_test/lib/python3.8/site-packages/shap/utils/_clustering.py:69: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
##   def _reverse_window_score_gain(masks, order, start, length):
## /home/hrodenhizer/miniconda3/envs/planet_processing_test/lib/python3.8/site-packages/shap/utils/_clustering.py:77: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
##   def _mask_delta_score(m1, m2):
## /home/hrodenhizer/miniconda3/envs/planet_processing_test/lib/python3.8/site-packages/shap/links.py:5: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
##   def identity(x):
## /home/hrodenhizer/miniconda3/envs/planet_processing_test/lib/python3.8/site-packages/shap/links.py:10: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
##   def _identity_inverse(x):
## /home/hrodenhizer/miniconda3/envs/planet_processing_test/lib/python3.8/site-packages/shap/links.py:15: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
##   def logit(x):
## /home/hrodenhizer/miniconda3/envs/planet_processing_test/lib/python3.8/site-packages/shap/links.py:20: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
##   def _logit_inverse(x):
## /home/hrodenhizer/miniconda3/envs/planet_processing_test/lib/python3.8/site-packages/shap/utils/_masked_model.py:363: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
##   def _build_fixed_single_output(averaged_outs, last_outs, outputs, batch_positions, varying_rows, num_varying_rows, link, linearizing_weights):
## /home/hrodenhizer/miniconda3/envs/planet_processing_test/lib/python3.8/site-packages/shap/utils/_masked_model.py:385: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
##   def _build_fixed_multi_output(averaged_outs, last_outs, outputs, batch_positions, varying_rows, num_varying_rows, link, linearizing_weights):
## /home/hrodenhizer/miniconda3/envs/planet_processing_test/lib/python3.8/site-packages/shap/utils/_masked_model.py:428: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
##   def _init_masks(cluster_matrix, M, indices_row_pos, indptr):
## /home/hrodenhizer/miniconda3/envs/planet_processing_test/lib/python3.8/site-packages/shap/utils/_masked_model.py:439: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
##   def _rec_fill_masks(cluster_matrix, indices_row_pos, indptr, indices, M, ind):
## /home/hrodenhizer/miniconda3/envs/planet_processing_test/lib/python3.8/site-packages/shap/maskers/_tabular.py:186: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
##   def _single_delta_mask(dind, masked_inputs, last_mask, data, x, noop_code):
## /home/hrodenhizer/miniconda3/envs/planet_processing_test/lib/python3.8/site-packages/shap/maskers/_tabular.py:197: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
##   def _delta_masking(masks, x, curr_delta_inds, varying_rows_out,
## /home/hrodenhizer/miniconda3/envs/planet_processing_test/lib/python3.8/site-packages/shap/maskers/_image.py:175: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
##   def _jit_build_partition_tree(xmin, xmax, ymin, ymax, zmin, zmax, total_ywidth, total_zwidth, M, clustering, q):
## /home/hrodenhizer/miniconda3/envs/planet_processing_test/lib/python3.8/site-packages/shap/explainers/_partition.py:676: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
##   def lower_credit(i, value, M, values, clustering):
## The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
## The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.

Prep Google Drive Authentication

Define Functions

assign_conf_stars

avg_precision

bbox

eval_expression

get_rast_id

googledrive_download

input_as_df

import_pred

lm_contrasts

lm_summary

make_my_dir

my_lm

plot_prediction

plot_zonal_stats

pred_as_poly

recall_precision

trim_outliers

val_as_poly

Prep Plot Variables

Load Data

Polygons

## Reading layer `rts_polygons_for_Yili_May_2022_v2' from data source 
##   `/home/hrodenhizer/Documents/permafrost_pathways/rts_mapping/rts_data_comparison/data/rts_polygons/rts_polygons_for_Yili_May_2022_v2.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 138 features and 11 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -139.218 ymin: 68.99102 xmax: 124.304 ymax: 72.99794
## Geodetic CRS:  WGS 84

Download Files

## [1] "There are 69 prediction tiles."
## [1] "There are 138 input data tiles."

Maxar GeoTiffs

Planet GeoTiffs

Sentinel GeoTiffs

Get Tile Bounding Boxes

Convert Predictions to Vector

Join Prediction Polygons into polys SF Dataframe

Convert Validation to Vector

Join Validation Polygons into polys SF Dataframe

Convert Input Data to DF

Interactive Map of Features

IoU

Calculate Intersection and Union

Calculate IoU

IoU (Calculated in Model) |Imagery|Test IoU|Val IoU| |:—:|:—:|:—:| |Maxar|0.73|0.72| |Planet|0.71|0.70| |Sentinel-2|0.68|0.66|

IoU (Calculated Directly from Predictions)
Imagery Val IoU
Maxar 0.67
Planet 0.64
Sentinel-2 0.61

Mean Average Precision

Recall-Precision Curves

MAP

Plot

Precision measures false positives. Recall measures false negatives.

This figure will only be interesting once we have added in negative training data.

Performance by Feature Size

Calculate Area

Size Distribution by Region

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

yg mean_size min_size max_size median_size
Other 20496.33 484 107280 10380
Yamal/Gydan 11409.41 512 47696 6204

Plot

Raw IoU Scores:

This is complicated by the fact that the rts_area column is calculated from the raster validation layer, which may contain several RTS features within one tile. Use rts_area (from the original RTS delineation), instead.

Run nls models and bootstrap parameters

Bootstrap predictions for plotting the nls models

Plot the Size/Performance plot

Active vs. General RTS Performance

imagery p_val x_pos star_y_pos label_y_pos p_label star_label
Maxar 0.1325867 1.5 0.9 0.95 p-value = 0.133
Planet 0.7955779 1.5 0.9 0.95 p-value = 0.796
Sentinel-2 0.5639142 1.5 0.9 0.95 p-value = 0.564

Plot

It is possible to get rid of the inner panel borders, if I decide that looks better: https://stackoverflow.com/questions/46220242/ggplot2-outside-panel-border-when-using-facet

Drivers of Unexpected RTS Prediction Performance

Classify Features Using Confidence Interval Approach

This approach first uses the 95% CI of the model parameters to determine whether RTS features were predicted better or worse than expected based on the model. Next, the threshold at which RTS size doesn’t impact IoU is determined from where the slope of the model approaches 0 (currently using slope < 1e-06) for each imagery type. RTS features smaller than this threshold that were predicted better than expected are analyzed later to determine why some small RTS can be identified from the imagery.

Is performance class the same across imagery sources for most RTS features?

Zonal Statistics

These plots summarize the input data values (mean or standard deviation) in RTS cells, background cells, and the normalized difference between the two (Delta = RTS - Background). Most of the input layer names should be self explanatory, but for the others:

lum = luminance = 0.299*r + 0.587*g + 0.114*b

sr = shaded relief

A few takeaway points:

  • NIR and NDVI (mean) figures indicate that RTS are predicted poorly where background cells have low NDVI/plant growth (plants reflect NIR, so have a higher NIR value). I.e. it’s harder for the model to find the RTS where there’s not much plant growth across the landscape.

  • NDWI (mean) figure indicates that RTS features are predicted better where it’s drier (fewer lakes?).

  • The standard deviation of elevation figure indicates that where there is higher variability in elevation within RTS relative to the background cells, RTS features are identified more easily. I.e. smooth landscape with uneven RTS is easy for the model to find.

I didn’t find much interesting about the rest of these figures, but they are included here for completeness:

Shapley Additive Explanation Analysis

Read in Models from Google Drive

## Automated hyper-parameter determination is applied with the following details:
## ----------
##  Number of convolution filters after each full-scale skip connection: filter_num_skip = [32, 32, 32, 32]
##  Number of channels of full-scale aggregated feature maps: filter_num_aggregate = 160
## ----------
## deep_supervision = True
## names of output tensors are listed as follows ("sup0" is the shallowest supervision layer;
## "final" is the final output layer):
## 
##  unet3plus_output_sup0_activation
##  unet3plus_output_sup1_activation
##  unet3plus_output_sup2_activation
##  unet3plus_output_sup3_activation
##  unet3plus_output_final_activation
## <tensorflow.python.checkpoint.checkpoint.CheckpointLoadStatus object at 0x7f393010bb50>
## Automated hyper-parameter determination is applied with the following details:
## ----------
##  Number of convolution filters after each full-scale skip connection: filter_num_skip = [32, 32, 32, 32]
##  Number of channels of full-scale aggregated feature maps: filter_num_aggregate = 160
## ----------
## deep_supervision = True
## names of output tensors are listed as follows ("sup0" is the shallowest supervision layer;
## "final" is the final output layer):
## 
##  unet3plus_output_sup0_activation
##  unet3plus_output_sup1_activation
##  unet3plus_output_sup2_activation
##  unet3plus_output_sup3_activation
##  unet3plus_output_final_activation
## <tensorflow.python.checkpoint.checkpoint.CheckpointLoadStatus object at 0x7f3899d79f70>
## Automated hyper-parameter determination is applied with the following details:
## ----------
##  Number of convolution filters after each full-scale skip connection: filter_num_skip = [32, 32, 32, 32]
##  Number of channels of full-scale aggregated feature maps: filter_num_aggregate = 160
## ----------
## deep_supervision = True
## names of output tensors are listed as follows ("sup0" is the shallowest supervision layer;
## "final" is the final output layer):
## 
##  unet3plus_output_sup0_activation
##  unet3plus_output_sup1_activation
##  unet3plus_output_sup2_activation
##  unet3plus_output_sup3_activation
##  unet3plus_output_final_activation
## <tensorflow.python.checkpoint.checkpoint.CheckpointLoadStatus object at 0x7f38993e9220>

Access model input data from Python